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ABSTRACT 

The internal angular momentum distribution of a star is key to determine 
its evolution. Fortunately, the stellar internal rotation can be probed through 
studies of rotationally-split nonradial oscillation modes. In particular, detection 
of nonradial gravity modes (g modes) in massive young stars has become feasible 
recently thanks to the Kepler space mission. Our aim is to derive the internal 
rotation profile of the Kepler B8V star KIC 10526294 through asteroseismology. 

We interpret the observed rotational splittings of its dipole g modes using four 
different approaches based on the best seismic models of the star and their rota¬ 
tional kernels. We show that these kernels can resolve differential rotation within 
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the radiative envelope if a smooth rotational profile is assumed and the observa¬ 
tional errors are small. Based on Kepler data, we find that the rotation rate near 
the core-envelope boundary is well constrained to 163 ±89 nHz. The seismic data 
are consistent with rigid rotation but a profile with counter-rotation within the 
envelope has a statistical advantage over constant rotation. Our study should 
be repeated for other massive stars with a variety of stellar parameters in order 
to deduce the physical conditions that determine the internal rotation profile of 
young massive stars, with the aim to improve the input physics of their models. 

Subject headings: Asteroseismology - Stars: rotation - Stars: oscillations (in¬ 
cluding pulsations) - Stars: individual: KIC 10526294 


1. Introduction 


One of the major important ingredients in the computation of stellar evolution models, 
from the star formation process until the death of the star, is rotation (e.g., 
for a recent monograph in the subject). Even in the case of slow rotation, the dynamical and 
mixing processes related to it are not negligible as they substantially affect the structure of 
the star. Unfortunately, direct measurements of the internal rotation profile of stars are not 
possible. Hence the inclusion of rotational effects in models rests on uncalibrated theoretical 
prescriptions. 


Maeder 2009 


It is a fortunate circumstance that the oscillation frequencies of a star are affected by its 


rotational properties (Ledoux 1951). The frequency splitting of the oscillation modes is well 


understood in the case of slowly rotating stars for which a first-order perturbation method 


is sufficient to model the frequencies (e.g. Aerts et al. 2010, for an extensive description). 
In this work, we assume that we are dealing with an unevolved star that does not possess 
a magnetic field and whose central frequencies of the rotationally-split multiplets are not 
affected by the slow rotation. Moreover, we assume that the deformation from spherical 
symmetry due to the centrifugal forces can be ignored. In that case the frequency splitting 
of the oscillation modes, as measured in the observer’s frame, is due to a combination of 
mode advection and the Coriolis force and can be computed from the so-called rotational 


kernels (cf. Eq.(3.356) in Aerts et al. 2010) 


Helioseismology delivered a very detailed view of the internal rotation profile fl(r, 6) of 
the Sun for the radial range r e [0.2 R Q , 1.0 R & ] ( R 0 denotes the solar radius) and for all co¬ 
latitudes 9, through frequency inversion of the rotational splittings of its hundreds of detected 


acoustic modes (e.g., Christensen-Dalsgaard 2002; Thompson et al. 2003), a technique that 
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has found its way even to laboratory experiments (Triana et al. 2014). Given that acoustic 
modes do not have sufficient probing power in the very inner regions and that the Sun does 
not reveal gravity modes, it is not possible to deduce the rotational profile for r < 0.2 R Q . 


ft is not currently within reach to derive the rotation profile for distant stars with 
similar precision as the Sun, but applications of asteroseismology did allow to deduce aver¬ 
aged rotation rate ratios D core /fl enV eiope from forward modeling of the rotational splitting for 


a few main-sequence B stars from ground-based monitoring campaigns (Aerts et al. 2003 


Pamyatnykh et al. 

2004 

Briquet et al. 

2007 

), as well as two 5 Sct-q Dor-type hybrid main- 


Moreover, f2 core was derived from gravity-dominated mixed modes in hundreds of red giants 


observed with Kepler (Beck et al. 2012; Mosser et al. 2012), while the estimate of their 


^envelope is uncertain due to the remaining dominant influence of the core regions on the 
measured splittings of the pressure-dominated mixed modes. Deheuvels et al. (2012) and 


Deheuvels et al. (2014) performed frequency inversions for seven selected subgiants in differ¬ 


ent evolutionary stages relying on the splitting of their dipole mixed modes. They selected 
profiles representing a linear decrease in rotation frequency in the core regions followed by 
a constant rotation profile in the extended convective envelope, with f2 core /flenveiope ranging 
from 2 to about 20. A similar result was obtained for the red giant KIC 5006817, which is 
the primary of an eccentric binary, by Beck et al. (2014). Only in some of those studies, 
e.g. Deheuvels et al. (2014), statistical model comparison was used to evaluate the likelihood 


of the optimal shape of the rotational frequency throughout the stars and few continuous 
and discontinuous functions for fl(r) were considered. Despite their frequent occurrence in 
e.g., geophysics, counter-rotating solutions were considered inappropriate for stars. Indepen¬ 
dently of that restriction, the Kepler results for D(r) so far deliver an important calibration 
for the improvement of evolutionary models for single and binary low-mass stars, given that 
the theoretical predictions of angular-momentum transport result in core rotation rates are 


higher by at least an order of magnitude than observed (see e.g. Eggenberger et al. 2012 


van Saders & Pinsonneault 2013 Cantiello et al. 2014 for the input physics in question). 


The need for redistribution and loss of angular momentum during evolution is also required 
on the basis of the internal rotation properties of white dwarfs derived from both forward 


modeling and inversion of their rotationally-split g mode oscillation frequencies (Charpinet 


et al. 2009; Corsico et al. 2012). 


In this work, we provide the first frequency inverted rotation profile of an unevolved 
intermediate-mass B-type main-sequence star from its rotationally-split dipole gravity modes 
detected in four years of Kepler data. The paper is organized as follows: we summarize the 
observational data and the resuls of forward seismic modeling in Section 2. In Sections 3 
and 4 we examine in detail the rotational kernels and the splittings associated with linear 
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rotation models. In Section 5 we present piece-wise two- and three-zone rotation models. 
Section 6 is devoted to inversion methods, including theory and results, and in Section 7 we 
present results for the rotation profile based on Monte Carlo simulations. We summarize 
and conclude in Section 8. 


2. Observational input and results of forward seismic modeling 


A first detailed asteroseismic analysis of the main sequence B star KIC 10526294, was 


presented by Papics et ah (2014), including an estimation of the amount of core overshooting 


following earlier approaches for main-sequence B stars with a well-developed convective core. 


They characterized KIC 10526294 as a slowly rotating SPB star (see e.g. Aerts et al. (2010) 
for a definition) exhibiting a series of 19 quasi-equally spaced dipole modes. KIC 10526294 
is so far the only multiperiodic SPB star with unambiguous detection of rotationally-split 
triplets from the Kepler light curve. For this reason, it allows us to probe its interior structure 
to a deeper level than for any other SPB so far. 


With the purpose of detailed seismic modeling, Moravveji et ah (2015) elaborated on 


optimal frequency error estimation, taking into account the signal-to-noise ratio, sampling, 


and correlated nature of the Kepler data, following the method by Degroote et ah (2009). 


This resulted in a correction factor of 3.0 to be applied to the formal errors obtained from the 
nonlinear least-squares fit. We estimated the splitting for each dipole mode as the average 
splitting between the measured m — +1 and m — — 1 peaks with respect to the central m — 0 
peak. This comes down to considering only the symmetric component of the splittings. The 
total variance was then estimated as the variance of the symmetric component plus the 
znter-variance: 

( 1 ) 


2 

°Total 


— - (cr\ + <y\ i) + cr 2 (5_i, h+i), 


where h-ti denotes the m = ±1 splittings and a±\ their individual uncertainties. This results 
in larger errors for those splittings with larger asymmetric components. Further, as explained 
by the authors, the error estimates for the triplet components as computed from the Rayleigh 


limit and taking into account mode crowding effects are too large (Papics et al. 2014, their 


Fig. 8). Nevertheless, we also used those overestimated values with the argument that they 
deliver the most conservative upper limit to the true frequency errors. We list both error sets 
in Tableland show in this work that our conclusions on ff(r) are essentially independent on 
the choice of error set. The best error estimates (Error Set 1) are used throughout the main 
text, while all results relying on the too large errors (Error Set 2) are treated in Appendix|A] 


Following an essentially identical approach as in |Papics et al. (2014), Moravveji et al. 


(2015) were able to End seismic models that match more closely the observed frequencies 
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Table 1. Symmetric components of the observed rotational splittings of KIC 10526294 


Central frequency [/iHz] 

Splitting 6 n im [nHz] 

Error Set 1 [nHz] 

Error Set 2 [nHz] 

5.4655 

45.28 

14.72 

16.54 

5.6272 

29.49 

0.69 

8.70 

5.7978 

33.91 

7.37 

8.81 

5.9873 

32.64 

7.47 

10.33 

6.1739 

41.99 

4.42 

7.56 

6.3959 

35.74 

0.36 

4.91 

6.6200 

29.43 

4.06 

7.38 

6.8703 

30.42 

0.55 

8.55 

7.1235 

33.07 

3.46 

6.09 

7.4213 

29.99 

1.00 

9.41 

7.7616 

41.55 

15.04 

17.17 

8.1163 

28.73 

1.03 

6.42 

8.5036 

29.50 

2.59 

7.75 

8.9398 

28.18 

3.18 

7.23 

9.4090 

27.53 

0.85 

4.63 

9.9115 

26.11 

1.67 

5.82 

10.4495 

26.41 

7.02 

10.04 

11.0429 

25.74 

0.49 

5.21 

11.7293 

23.32 

4.55 

9.34 
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Fig. 1.— Top: Observed rotational splittings (top panel, blue error bars), and the symmetric 
component used as inversion inputs (top panel, red error bars). The errors on the inversion 


inputs were taken from Moravveji et al. (2015) and the green error bars are from Papics 


et al. (2014), see text for details). We chose the abscissae of the inversion inputs so as to 


coincide with the central m = 0 peak of the observed triplets. Bottom: the location of the 
central peaks as black dots, and the mode periods from models: vertical solid lines for the 
best model in Moravveji et al. (2015) (Model 1), vertical dashed lines for the best model 


m 


Papics et al. (2014) (Model2). 


thanks to the inclusion of extra diffusive mixing in the stellar envelope, in addition to core 
overshooting. Such additional mixing was already found necessary for the B3V SPB star 
HD50230 (Degroote et al. 2010). The parameters of the two models are given in Tablc[2j 
We base the present work on the best matching Model 1 from Moravveji et al. (2015) but we 
checked that the main qualitative features of the resulting rotation profiles from inversion do 
not depend on the choice of the seismic model, as long as it is able to reproduce reasonably 
well the main observed characteristics of the star. We return to this in Section liOl 


In a non-rotating star, oscillation mode families are characterized by their radial order 
n and degree /, with individual family members corresponding to different values of m, the 
azimuthal wave number, sharing the same eigenfrequency. Rotation lifts this degeneracy. 
The identification of the radial order n of the 19 detected modes of KIC 10526294 was 
achieved by comparing the periods of the observed zonal (m = 0) dipole (Z = 1) g modes 
with those predicted by equilibrium models computed with the MESA stellar structure and 
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Fig. 2.— The rotational kernels of dipole zonal modes of highest and lowest radial order, 
plotted against the stellar fractional radius r/R*, for the best matching models found by 
Papics et ah (2014, lower panel) and Moravveji et ah (2015, upper panel) — cf. lower panel 


of Fig.[lJ Some of the modes in the upper panel are trapped near the convective core as 
evidenced by the large peak just outside the convective core, as visible in the inset. The 
kernels from the best seismic model found by Papics et al. (2014) shown in the bottom panel 
do not exhibit trapping, but the overall shapes are similar otherwise. Inversion results from 
both of these models are qualitatively similar. 


evolution code ( 

Paxton et al. 

2011 

2013 

) and coupled to the GYRE stellar oscillation code 

(Townsend & Teitler 

2013 

). This comparison is possible thanks to detected triplets on the 


one hand, and the almost equally-spaced sequence of dipole g modes with consecutive order 
as theoretically expected for such modes on the other hand. The radial orders matching the 
observations range from n = 32 to n — 14. 


Most of the 19 dipole modes of KIC 10526294 reveal a very narrow rotationally-split 
triplet structure. If we assume that the cyclic rotation frequency, which we denote as 12, 
depends on the radial coordinate r only, the frequency splitting of a mode with degree l and 
radial order n, denoted as S n i m , can be written as 

r R * 

8nim = m/3 n i K ni (r) 12(r) dr, (2) 

Jo 

where the unimodular mode kernel K n i(r) is a function of the mode’s displacement ampli- 



































































tudes £ r (r) (vertical) and Ch( r ) (horizontal), while /3 n i is given by 


y l = C [£r + 


2£rih - Q r 2 pdr 


C [Cr+K l + 1 )Qr 2 pdr 


and is connected with the Ledoux splitting as j3 n i — 1 — C n i (e.g., Aerts et al. 2010, Chapter 
3). The rotationally-induced splittings for the 19 detected modes are shown in Fig. 0 (top 
panel, adapted from Fig. 8 in Papics et al. (2014)), together with the mode periods derived 


from Model 1 and Model2 (bottom panel). The mode kernels for both models are shown in 
Fig .[2] and their squared Brunt-Vaisala frequency N 2 is shown in Fig.|3} The latter illustrates 
the slightly more advanced stage of Model 1 compared to Model 2. 

Some of the observed splittings are not symmetric with respect to the central m — 0 
peak. These asymmetries are usually related to mechanisms capable of lifting partially 
the 21 + 1 degeneracy of a given multiplet, such as deviations from sphericity or large- 
scale magnetic fields. We assume to be dealing with a physical phenomenon that causes 
asymmetries without lifting the degeneracy between the retrograde (m < 0) and prograde 
(m > 0) modes. Rotationally-induced splittings as expressed by Eq. ([ 2 ]) are only capable of 
lifting the ±m degeneracy. Hence, forward modeling including only differential rotation up 
to first order cannot fully capture these observed asymmetries. On the other hand, Eq. (J2| is 
still perfectly valid even in the presence of mechanisms mentioned above, with the caveat that 
it then accounts for one-half of the splitting between the +m and the — m modes and not for 
the splitting between a mode with a given m and the corresponding central m = 0 peak. It is 


known that the presence of a magnetic field near the convective core, as discussed by Hasan 


et al. (2005), can give rise to such effects. In fact, the mode kernels considered here have 


substantial amplitudes precisely in that zone, making them particularly susceptible to this 
effect. Given that we have no information on the presence or absence of a magnetic field in 
KIC 10526294, we will address only the symmetric components of the splittings in this study, 
assuming rotation to be the dominant mechanism responsible for the observed splittings and 


Table 2. Fundamental stellar parameters of KIC 10526294 from the best matching 


theoretical model (Model 1) found by Moraweji et al. (2015) and by Papics et al. (2014) 


(Model 2). 


Model 

Teff [K] 

M*/M 0 


fov 

Z 

W 

Age [Myr] 

x 2 

1 

13000 

3.25 

2.215 

0.017 

0.014 

0.627 

63 

1.42 

2 

12470 

3.20 

2.100 

< 0.015 

0.020 

0.693 

12 

10.9 
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leaving the modeling of the asymmetries for a future more specialised study. As explained 
earlier, the presence of the asymmetries leads to increased errors in the splittings. 



Fig. 3.— The squared Brunt-Vaisala frequency N 2 as a function of stellar fractional radius 
for the best model from Moravveji et al. (2015) (continuous blue curve) and for the best 
model from Papics et al. (2014) (red-dashed curve). The larger peak of the former model 
near the core boundary (at r ~ 0.15 i?*) is comparatively further away from that boundary 
than the peak of the latter model, thus allowing some modes to be “trapped” between the 
core boundary and the peak. See also Fig. [2} 


3. Cumulative kernel integrals 

Each mode samples differently the internal rotation of the star. This is usually illustrated 
by plotting the cumulative integral of the kernels JFj. For increased contrast we plot instead 
in Fig. [4] the cumulative integral of fcj(r), which we define by 

ki(r) = Ki{r) - (K{r)) , (3) 

where ( K(r )) is the average of the kernels across modes at each radius r. This (K{r )) is 
the ‘common’ kernel and its integral will only contribute to the average of the splittings. 
Similarly we can express a given profile f2(r) as a sum of its mean value (across the radial 
coordinate), plus a fluctuating part uj(r) (with zero mean): 

f2(r) = Cl + c o(r). (4) 

Therefore we can write the scaled splittings as: 

A i = nl ™ =Q+ (K(r)) (j(r) d r+ ki(r)cj(r) dr. (5) 

m Pnl JO Jo 
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The first two terms on the right hand side comprise the average splitting, while differences 



Fig. 4.— Cumulative integrals of ki{r ) (top) and r ki(r) (bottom) for a sample of dipole 
(/ = 1) modes from Model 1 with various degrees of trapping. Mode with n — 27 is the 
“most” trapped near the core as compared with the n — 21 which is the “least” trapped. 

across modes come into play in the last term. Figure[4] shows the cumulative integrals of 
fci(r) and of fcj(r) r for a sample of modes from Model 1 with various degrees of trapping near 
the core. From the figure it is evident that the “least” trapped mode (the one with n = 21) 
would have the largest splitting in response to a rotation profile that increases linearly with 
radius. Conversely, the “most” trapped mode near the core (n = 27) would have the smallest 
splitting under the same condition. 
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4. Trapped modes and linear rotation profiles 


Here we follow closely the analysis done by Kawaler et al. (1999) for g inodes in white 
dwarf pnlsators. We take advantage of the fact that some modes are trapped as revealed 
by the kernel amplitudes. They are trapped very close to the overshooting zone while other 
modes have more spread-out amplitudes, comparatively. This trapping manifests itself as 
reduced period spacings if we plot them as functions of period (Kawaler & Bradley 1994). 
The observed spacings of KIC 10526294 are shown in Fig.[5]together with the scaled splittings 
Kim./m(3 nl (see Eq. §• A linear fit to the scaled splittings results in a slope of 27.04 nHz/(10 5 s) 
and an intercept of 26.20 nHz. 



Period [10 5 s] 


Fig. 5.— The scaled rotational splittings 5 n i m /mf3 n i deduced from the observations (red 
error bars) as a function of the mode period together with a linear fit (dotted red line, slope: 
27.04nHz/(10 5 s), intercept: 26.20nHz, scale is on the left). The dashed line connects the 
observed period spacings defined as A P n = n n — n n _!. The black continuous line connects 
the spacings derived from the seismic Model 1 in Table 1 while the grey continuous line 
connects the spacings derived from Model 2. (scale is on the right). The trapped modes of 
Model 1 are those with periods close to 0.95 x 10 5 s and 1.6 x 10 5 s. 

Let us now do forward modeling to find the predicted splittings using synthetic, linear 
rotation test profiles. Two test profiles have H = 5 nHz at r — 0 and have the same slope in 
absolute value, 1 nHz /R*, but opposite in sign. The results are shown in Fig.[6j where the 
top and bottom plots correspond to increasing and decreasing rotation profiles, respectively. 
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We see the clear signature of mode trapping. Trapped modes are closer to the core, so if the 
rotation profile increases with radius, then the corresponding splittings are comparatively 
smaller. Analogously, if the rotation decreases with radius then the trapped modes will 
show comparatively larger splittings than the other modes. The latter situation is precisely 
what we can see in Fig. [6} A similar situation occurs for the white dwarf PG 1159-035, as 
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Fig. 6.— Top: the scaled rotational splittings for a linearly increasing rotation profile derived 
from Model 1 and 2 (continuous blue and grey curves, respectively), the blue dotted line is 
a linear fit to the splittings from Model 1, scale is on the left. Bottom: the splittings for a 
linearly decreasing rotation profile derived from both Model 1 and 2 (continuous green and 
grey curves, respectively), the green dotted line is a linear fit to the splittings from Model 1, 
scale is on the left. For comparison, the horizontal lines show the splittings from constant 
rotation profiles having the same mean splittings as the splittings from the increasing or 
decreasing profiles (top, blue horizontal line and bottom, green horizontal line, respectively). 
The black dashed lines on both plots show the period spacings from observation, scales are 
on the right. 


reported by Kawaler et al. (1999), the only difference being that some modes in the white 
dwarf are trapped close to the surface such that the results are reversed compared to those 


for KIC 10526294. 


This simplified analysis is helpful because it gives a first idea of the sign of the slope of 
an unknown rotation profile just by plotting the splitings and the period spacings and see if 
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they vary in phase. In our case it is not very clear if they vary in phase or not, so instead we 
make use of the linear trends of the splittings as indicated with the linear fits (dotted lines) 
in Fig. [6} When period spacings and splittings vary in phase, the linear trend is downward 
(negative slope). Conversely, when period spacings and splittings are in anti-phase, the 
linear trend is upward (positive slope). The observed splittings of KIC 10526294 have an 
increasing trend, as evidenced by the linear fit (dotted red line) in Fig. [5] We associate this 
with a decreasing rotation rate. 

We can now do a simple calculation to estimate the optimal slope a of a linear rotation 
profile Q(r) = a r + b that best matches the observed slope. We assume that the slope of the 
fit to the measured splittings is linear with respect to the slope of the rotation profile, which 
seems to be the case. The ratio of the linear slopes associated with the observed splittings 
and the splittings of the linearly decreasing test profile discussed above is ~ 2.0918 x 10 3 . 
Therefore, if we use a linearly decreasing rotation profile with a slope of 2.0918 x 10 3 times 
the original slope of our (linearly decreasing) test profile and adjust its mean level to match 
the mean of the observed splittings, we might get an idea of the underlying rotation profile. 
The slope we obtain in this way is a ~ — 2.0918/iHz /R*. Adjusting the slope of the test 
profile so as to match the trend of the observed splittings also requires adjusting b according 
to Eq. ([6]) if we are to match the average of the observed splittings. We calculate the intercept 
b as (see also Eq.[5]) 

rR* 

b = (6) — a / (K(r)) rdr, (6) 

Jo 

where () denotes the average over all modes, and get b ~ 0.957 yuHz. These values of a and b 
imply that part of the rotation profile becomes negative. Hence, in the case that the splittings 
are actually caused by a linear rotation profile or any other profile closely resembling a linear 
one, we deduce that the mean observed splittings should have been considerably higher in 
magnitude in order to obtain a rotation profile that would not include counter-rotation inside 
the star. Of course, with this exercise, we were only trying to match the linear trend and the 
mean of the observed splittings through a linear test profile. In this case the rms deviation 
of the predicted splittings from the observed ones is around 22.42 nHz, comparatively large 
compared to the mean error from the observations which is 8.41 nHz (scaled, from Error 
Set 1). We can make use as well of the reduced x 2 values which we compute throughout this 
work as 



where Si and S t are the measured and predicted splittings, respectively, e* are the errors, u is 
the effective number of degrees of freedom and M is the number of observed splittings. In 
the case just discussed above we fitted two parameters, a and b, so the number of degrees of 
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freedom is v = M — 2 = 17 which leads to y 2 ~ 194.8. 

Figure[7] presents another hint towards the decreasing rotation as the radins increases. 
We plot the observed splittings as a function of k t (r) r dr, where ki(r ) is defined in 
Eq.g. This should have been a straight line with negative slope if the splittings were 
actually caused by a linear rotation profile. As the figure shows, the splittings have an 
overall downward trend which is a rough indication that the rotation profile decreases as we 
move towards the star’s surface. 
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Fig. 7.— The observed splittings (scaled) as a function of ki(r)rdr. Even considering 
the two splittings marked in red (which are coincidentally the ones with largest errors), there 
is a downward trend which hints at decreasing rotation as the radius increases. 


n=32 

Period = 1.83 x 10 5 s | 


n=22 

Period = 1.29 x 10 5 s I 


o° 


In a next exercise, we searched for the linear profile that minimizes y 2 . This profile has 
a slope of a = — 530.7 nHz/i?* and an intercept b = 285.1 nHz, leading to y 2 = 6.74. This 
profile again leads to negative values for the rotation frequency in the outer envelope. To have 
an idea of the comparative statistical significance of this result, we computed the optimal 
constant rotational profile, as well as the optimal linear profile restricted to positive values 
(using a Lagrange multiplier as an additional fitting parameter). We obtained y 2 onst = 19.86 
and x+ — 16.86, respectively. None of the positive linear rotational profiles have a y 2 similar 
to the linear profile with counter-rotation. 
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Period [10 s s] 


Fig. 8.— Top left: x 2 from the difference between the predicted and observed splittings as 
a function of the two-zone parameter tq. Top right: the values of Oi (blue) corresponding 
to r < ro and f2 2 (green) for r > r^, both as functions of r$. The vertical dotted line in 
black indicates the extent of the core and the overshoot region, it coincides closely with 
the optimum r 0 . Bottom left: y 2 computed from a three-zone model with rq and r 2 as 
parameters. The lowest x 2 is achieved when r\ = 0.85 i?* and r 2 = 0.91 i?*. Bottom right: 
Observed and predicted splittings from the two-zone model together with the predicted 
splittings corresponding to the lowest y 2 from the three-zone model (green squares). 

5. Linear, piece-wise rotation models 

We now assume a two-zone, piece-wise rotational profile such that its value is O] if 
0 < r < r 0 and 0 2 if r 0 < r < R*. The parameter r 0 is variable and we optimize x 2 ( r o) so 
as to best match the observations, whose errors are derived from Error Set 1. The resulting 
X 2 versus ro is shown in the top left panel of Fig.[8j The minimum occurs at vq ~ 0.166 i?*. 
The values of fli,f2 2 are shown as functions of ro in the top right panel of Fig.[8j The 
bottom right panel of Fig.[8] shows the observed splittings (scaled) and the splittings for the 
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two-zone model with the minimum at r 0 = 0.166 7?*, corresponding to Qi = 262.71 nHz and 
fi 2 = 49.33 nHz. We note that the boundary of the convective core of Model 1, extended 
with the core overshoot zone, is situated at 0.1652 7?* (Moravveji et al. 2015). The second 
deepest minimum in the top left panel of Fig. [8] occurs at r 0 = 0.741 7?* and corresponds to 
Hi = 142.0 nHz and H 2 = —418.2 nHz; this solution leads to an that does not play a 
special role in Model 1 in terms of physical quantities. 


The two-zone model thus favors a region rotating with a period of 42 d near the core 
overshoot zone and a co-rotating envelope with a period of 254 d (y 2 = 1.59). The second 
minimum has y 2 = 2.36 and corresponds to a counter-rotating profile. The averaging kernels 
(for a definition, see Section 6) associated with the best two-zone model are presented in 
Fig.[9]and reveal that the outer zone averaging kernel probes mostly the radiative envelope, 
while the inner zone averaging kernel exhibits a large maximum just before reaching the 
core-envelope boundary and rapid oscillations around zero away from it. 
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Fig. 9.— Top: averaging kernels J^(r',r) of the optimum two-zone model with a disconti¬ 
nuity at ro/-R* = 0.166. Bottom: zoom on the core-envelope boundary region. 

A slightly more complex version of this two-zone model is implemented by introducing a 
third middle zone where the rotation profile changes linearly from fix at the end of the inner 
zone to fl 2 at the start of the outer zone. In this case, we have two linear parameters fli, h2 2 
and two non-linear parameters ri,r 2 defining the zone boundaries. The lowest x 2 = 1.51 
for this three-zone model is achieved when ry = 0.85 77* and r 2 = 0.917?* (see bottom left 
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panel of Fig.[8]) and indicates counter-rotation (fix = 151.8nHz, — —1069.5nHz). The 
corresponding splittings are also shown in the bottom right panel of Fig. [HI 


Very similar results are obtained when the uncertainties are derived from Error Set 2 or 
when the most asymmetric modes are excluded as evidenced by Figs. [19] through [21} All the 
three-zone y 2 minima in the four cases, i.e., using Error Set 1 or 2 and with or without the 
most asymmetric splittings, correspond to cases with counter-rotation. According to their 
y 2 values, these three-zone models have a statistical significance comparable with the best 
two-zone model. However, the position of their discontinuities rq and r 2 have no obvious 
physical meaning. 


To estimate the performance of different models, with the aim to choose the best one in 
a statistical sense, one can only rely on likelihood ratios in the case of nested models, i.e., for 
models where all terms of a simpler model version also occur in a more complex version of 
the model (e.g., Hastie et al. 2009). We are not in such a situation here because we wish to 
compare linear, discontinuous linear multi-zone, and continous non-linear inversion profiles 
(the latter will be discussed in the next Section). In such a case of non-nested models, an 
adequate statistical measure for model selection is the Akaike Information Criterion (AIC), 
which assigns a score to a given model rewarding goodness-of-fit (e.g., as measured by y 2 ) but 
penalizing overfitting, thus discouraging the use of complex models with too many adjustable 
parameters (Burnham & Anderson 2002 Hastie et al. 2009). For our purposes, having only 
a relatively small number of measurements (M = 19), it is appropriate to use the corrected 
AIC, denoted AICc, which we define according to its common use in the literature (e.g., 


Hurvich & Tsai 1989): 


AICc = 2 k + x 2 + 


2 k{k + 1) 


( 8 ) 


M — k — 1 ’ 

where k is the number of parameters to fit and M is the number of observations. As 
advocated by Burnham & Anderson (2002), k should include the variance of the residuals as 
a parameter to be fitted, e.g. k = 3 for a linear regression. With this definition, k = M — v+l. 
The preferred model among a set of models is the one with the lowest AICc value, where we 
limit proper model comparison to the case M — k > 1. These AICc values are only intended 
for model inter-comparison and have no absolute meaning by themselves. Table[3] lists all 
the rotation models considered in this work together with some of their associated statistical 
measures, including their AICc’s. Similar tables based on Error Set 2 or by avoiding the 
most asymmetric splittings, are given in Appendix[Aj We can see from the AICc values in 
Table[3]that the two-zone piece-wise model outperforms the three-zone model. 

We end this Section by noting that the overall spectral line broadening of 18kms _1 
measured for KIC 10526294, which is the combination of rotational and pulsational broaden¬ 
ing (Papics et al. 2014), is compatible with all the 0 2 -values found from the minima listed in 
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Table 3. Comparison of rotation profiles for Model 1 (except for one entry, where we used 
Model 2). The profiles marked with + are enforced to be positive-definite. The effective 
number of degrees of freedom v for the inversions were computed as M — k, where k is 
trace of the ‘hat’ matrix GC T (see Section |b|, following the same approach as Deheuvels 


et ah (2014) and explained in Hastie et ah (2009). Error Set 1 is used here. 


Rotation Prohle 

rms error [nHz] 

V 

x 2 

AICc 

Constant 

11.63 

18.00 

19.86 

362.3 

Linear 

12.32 

17.00 

6.74 

122.1 

Linear-I- 

11.46 

16.00 

16.86 

280.7 

Two-zone 

10.88 

16.00 

1.59 

36.36 

Three-zone 

9.48 

15.00 

1.51 

37.35 

RLS, N = 8 

9.31 

14.43 

0.52 

24.52 

RLS, N = 14 

10.46 

14.32 

1.44 

38.18 

RLS, N = 8 (Model 2) 

11.35 

14.74 

6.42 

110.3 

RLS+, N = 8 

12.05 

10.43 

17.16 

222.2 

RLS+, N = 14 

11.09 

5.32 

3.07 

184.3 

SOLA, N = 8 

9.77 

12.64 

0.58 

33.59 

SOLA, N = 14 

8.19 

9.28 

0.45 

60.11 
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this Section and does not allow any discrimination among those solutions as it was the case 
for the subgiant studied by Deheuvels et al. (2012). In the following, we investigate rotation 
profiles obtained from inversion methods. 


6. Inversions 


We first introduce the basic concepts and terminology behind the inversion approaches 
we applied. We are interested in the approximate determination of Q(r) based on a set of 


observed rotational splittings <5 n z m , see, e.g., Gough (1985) for one of the earliest applications 


of this method in the solar case and Kawaler et al. (1999) and Deheuvels et al. (2012) for 


applications to white dwarfs and subgiants, respectively. This constitutes a linear problem 
and the approximate solution f2(r) can be written as 


M 


£i(r) = y^c,(r)A,, 


(9) 


i —1 


where {i} represents the collective index {nlm}, A* = 8 n im/ml3 n i are the scaled splittings, 
Cj(r) are the yet-unknown inversion coefficients , and M is the number of observed modes. It 
is convenient to express the approximate rotational profile 0(r) in terms of the true profile 
fl(r) by means of the averaging kernels Jf(r',r). They are related to the kernels Kfr) 
through J^(r',r) = YlfLi c i{ r ')Ki{ r ) and fulfill 

r-R* 


fl(r') = / J^(r', r) fl(r) dr. 


( 10 ) 


From the preceding relation it is clear that the averaging kernels (r', r) should be as much 
as possible localized around r\ ideally resembling a delta function <5(r',r). 

We consider now in what follows a radial grid (scaled with the stellar radius Rf) of 
N + 1 uniformly spaced points {r 0 ,r!,... , rjv} with r 0 = 0 and r N = 1, covering the full 
range of fractional radius. The goal of the inverse problem is to determine the N unknowns 
Clj, which represent the predicted angular velocity f2 at radius r such that 7\,_i < r < rj, 
where j = 1,..., N is the grid index. This discretization of on a radial grid allows us to 
write an expression for the corresponding predicted splittings Aj, based on Eq. (|2]) , as: 

N f r :j 

A, : = Gjjtlj , where Gjj = / iG(r)dr, (11) 

j =i ^ r 9- 1 

or simply A = Gfl in matrix form. Analogously, we express Eq. (|9]) 

M 

vtj = 'y ] CjjAji 


as 


i =1 


( 12 ) 
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with the inversion coefficients c t -j constituting the matrix C. It is instructive to put the 
relation above in terms of the discrete version of the true rotation profile Q(r). It is straight¬ 
forward to show that the matrix A, defined through A = C T G, accomplishes such task by 
fulfilling 

N 

( 13 ) 

k =1 

Ideally Ajk should resemble a Kronecker-delta 8jk indicating that the recovered profile at 
a given radius (specified by the grid index j) does not suffer from ‘leakage’ coming from 
other radial regions. The matrix A is thus the discrete equivalent of the averaging kernels 
JC(r',r). 


Note that the observed splittings are linearly related to the predicted splittings through 
the matrix GC T , which is known as the ‘hat’ matrix. The trace of this matrix is an estimate 
of the effective number of adjustable parameters (Hastie et al. 2009, p. 232). 


If the observational errors e; are uncorrelated, as we assume here, the variance of the 
recovered profiles can be estimated as 


a 


A) = 


M 

i= 1 


(14) 


The relation above accounts for the errors on the measurements only and its impact on the 
inversion results; it does not account for the errors inherent to the inversion process itself. 

Below we describe two different inversion techniques we used to obtain an approximation 
to the internal rotation of KIC 10526294, as well as quantitative estimates of the uncertainties 
originating from the measurement errors. 


6.1. Regularized least squares method 


A technique commonly used in hclio- and asteroseismology is the regularized least 


squares (RLS or Tikhonov) method (e.g., Craig & Brown 1986), which seeks to minimize 
the quantity 2F defined as 


(At- A A 

& — ^ - -ji -hRLS 




i= 1 


e? 


d 2 n 

dr 2 


dr, 


(15) 


where A $ are the observed splittings, e t the corresponding measurement error, A,; are the 
predicted splittings, and fi rls is a, free parameter (known as the regularization or smoothing 
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parameter) used to limit the norm of the second derivative on the predicted f2(r). Using 
our discrete radial grid described earlier, minimization means 82?/dClj = 0 for all j. This 
condition can be written more explicitly using Eq. © and Eq. (fl5|) 


as 


H t (G Cl - A) 


hRLS 


L L S 7 = 0, 


(16) 


where Hij = Gij/ef. L being the discrete second derivative operator, 5 the radial grid spacing, 
and eg is a scale factor (introduced for convenience) which we set equal to the squared mean 
of the errors A formal solution is 


n = (^h t g + l t l^) h t a. (17) 

Obviously, the inverse matrix on the right hand side of Eq. ( [Tt| ) might not exist and in 
practice we seek a solution in the least squares sense instead. 


6.2. The SOLA method 


The Subtractive Optimally Localized Averaging (SOLA) method (Pijpers & Thompson 


1994) consists in determining a linear combination of the inversion coefficients c;(r') such that 
the averaging kernels resemble as much as possible a target function T(r',r) while keeping 
the variance of the predicted profiles, cr 2 (f2), low. To implement this method we minimize 




[2(2 (r', r) — T(r', r)] 2 dr + 


hSOLA 


-0 


M 


X) C ?( r ') e i 


(18) 


2—1 


at each r', with the constraint J^* Jf?(r',r) dr = 1. In addition to the free parameter /tsola 
we can also adjust the shape of the target function T. The problem reduces to solving the 
linear set of M equations {i — 1,..., M and for each radial location r') 


M t-R, 

Y, w * c ^ r ') = / K i( r ) T ( r '> r ) dr > 


(19) 


where W ik = f Q R * Kffr) K k (r) dr + (/xsolaAq) S ik ef, together with the constraint c k( r ') = 
1. Using the discrete radial grid with N segments and choosing the target functions as 


T(r', r) 


N if r' G (r-j-i, rf) 
0 otherwise, 
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the problem to solve becomes 


M 

E 

k=1 L 


I(GG T ) J + ^4f= 


Ckj NGij, 


( 20 ) 


with the constraint Y2t=i c kj — 1- I 11 this case the target f un ction will approach a Dirac-5 as 
N increases. On the other hand, N can be low or moderate as long as fl(r) can be assumed 
not to vary appreciably over a radial segment of the grid. 


6.3. Profiles from inversion 


We present the internal rotation profiles obtained by the two methods described above. 
For both the RLS and SOLA methods, we scanned different resolutions ranging from N = 2 
to 14, each covering a wide range of smoothing parameters /irls and Usola , to examine the 
resulting inversion profiles. To choose the appropriate parameters N and fi is not an easy 
task and in practice it depends on the a-priori information we might have. For example, if 
we can assume that the rotation profile does not change appreciably over a radial distance A, 
then we can safely use a resolution such that N = 1/A. Alternatively, we can interpret the 
inversion result as providing information only on those components of the rotation profile 
that do not change appreciably over a distance 1/N. Generally the higher the change of the 
profile over a given length scale is, the lower the accuracy of the inversion. If we consider that 
the rotation profile can be represented as a superposition of profiles with increasing detail 
(e.g. like a Fourier expansion), then the inversion methods can provide useful information at 
least on those components that change slowly with radius. Below we discuss the results for 
the N values with the best statistical score in terms of the AICc. 


We show in Fig. 10 the inversion profiles with /irls — 10 -5 and Usola = 10~ 2 , which 
represent a good balance between error and resolution for Error Set 1. Similar figures for 
Error Set 2 and/or ignoring the asymmetric splittings are shown in Figs. 22 24]in Appendix|A] 
We see that both methods give qualitatively similar results for N = 8 and N = 14. The 
uncertainty of the SOLA method grows quickly as N is increased, as opposed to the RLS 
method. The averaging kernels, as can be judged from Figs. [26] and 27 in Appendix|Bj are 
better localized using the SOLA method (as expected, by design), although the uncertainties 
are somewhat larger compared to the RLS method. The kernels feature larger amplitudes 
near r ~ 0.2 R*, corresponding to small uncertainties at that location. Further out in 
radius, near r ~ 0.9 R*, the uncertainty is larger but the rotation rate is still constrained 
to be opposite in sign. Note that the kernels provide no information for r < 0.15 R* so the 
inversion results within this radial range have no meaning (cf. Appendix[B]). An appropriate 
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Fig. 10.— Inversion profiles resulting from the RLS method (top row) with Hrls — 10 -5 and 
the SOLA method (bottom row) with usola = lO -2 for two different resolutions N = 8,14 
(left to right). Vertical dashed lines indicate the approximate location of the convective 
core’s radius. Error Set 1 has been used. 


assessment of the inversion’s accuracy is provided by the averaging kernels, or their discrete 
counterpart represented by the matrix A (see Figs. 26 and 27). Generally speaking, as [irls 
or I^sola increases, the predicted splittings A* will deviate more from the measured ones. 
The predicted splittings for N = 14 using both RLS and SOLA methods and Error Set 1 is 
displayed in Fig.[lTJ The result in the case of the omission of the asymmetric splittings is 
provided in Fig. [25] of Appendix |A[ 


The best two-zone model from Section 5 might seem in disagreement with the profiles 
obtained from inversion but upon further inspection they are actually in good agreement, at 
least for the outer zone. Indeed, if we take the profile from RLS inversion with N = 8, Hrls — 
10~ 5 and interpolate it appropriately, the resulting average rotation using the outer zone 
averaging kernel displayed in Fig. [9] amounts to 50.34 nHz, which is in very good agreement 
with the 49.33 nHz value obtained in Section 5 for the outer zone. The rotation rate computed 
from the inversion profile and the inner zone averaging kernel amounts to 206.01 nHz, which 
is to be compared with the 262.71 nHz value for the inner zone in the two-zone model. 
Although both the best RLS model and the best two-zone are mutually consistent, the RLS 
model resolves the outer zone better. The AICc values also give preference to the RLS 
inversions as we deduce from Tabled 
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Fig. 11.— The predicted splittings for both the RLS (top) and SOLA methods (bottom) 
using a resolution of IV = 14 and a range of /j parameters for each method. The predicted 
splittings were computed based on the uncertainties from Error Set 1. 


To add yet another model comparison, we have performed the so-called leave-one-out 
cross-validation technique (see e.g. Hastie et al. 2009), which does not rely on y 2 values. This 
technique consists in omitting one of the measurements when fitting a model and comparing 
the predicted splitting based on the fitted model A* with the measurement that has been 
omitted A f e* = (Aj — A^/e*. We performed this procedure for each of the measurements 
at a time and obtained a final score by computing the rms value of all the e*. The preferred 
model is then the one with the lowest score. In the case of the two-zone model this score is 
1.51 while for the best RLS inversion the score is 0.79. 


Using a synthetic profile without counter-rotation as a test, we found that the RLS 
inversion method works properly and that it is not prone to give spurious counter-rotation 
solutions. See Appendix [D] for further details. 
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It is possible to perform RLS regularized inversions enforcing a-priori a rotation profile 
that does not exhibit counter-rotation. We can achieve this by including additional terms to 
the quantity ST defined in Eq. (15) in such a way that they represent our prior knowledge of 
the rotation profiles being definite positive. Loosely speaking, we can express the probability 
of obtaining a value Qj at a given radial location as being proportional to e XjClj , with A j > 0. 
A concrete technique to solve such minimization problem is provided by the Karush-Kuhn- 


Tucker conditions (Kuhn & Tucker 1951), which consists in minimizing 5T—JT A jQj together 
with the constraints A jQj = 0 for all radial locations j (see e.g. Boyd &; Vandenberghe||2004). 


The A j are treated as Lagrange multipliers. The result of this exercise is shown in Fig. 12 



Fig. 12.— Left: comparison of inversion profiles obtained by restricting the rotation to be 
positive for all r and those using unrestricted RLS regularization for two different resolutions 
N = 8,14 and for Hrls = 1CT 5 . Right: the corresponding predicted splittings compared to 
observations. See main text for details. 


A qualitative idea of how well these restricted-positive profiles represent the data is 
provided by the plot on the right of Fig. 12 Quantitatively, the y 2 for the restricted-positive 
profiles are y| + = 112.9 and y 2 4+ = 46.1 corresponding to resolutions of iV = 8 and N = 14, 
respectively. These values are to be compared with y 2 = 0.52 and y 44 = 1.44 from the 
unrestricted case. Note here that the effective number of degrees of freedom v is reduced 
considerably when restricting the profiles to be positive definite given the additional param¬ 
eters included as Lagrange multipliers. In Table[3]we assembled all the AICc scores of these 
and other inversions along with results from previous sections. 


Although it is possible to increase the resolution N beyond the number of observations 
M without overfitting for the inversion result itself, following the principle of regularization 
(see Appendix JcJ , the value of v computed as the trace of the ‘hat’ matrix GC T might 
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become smaller than unity if the regularization parameters are small, pushing the AICc 
used for model comparison beyond meaningful values. For this reason we limit ourselves to 
moderate resolutions for the sake of meaningful model comparisons. The best models are 
the RLS or SOLA inversions with N = 8 with large statistical margin over all the other 
rotational profiles we obtained. 


The profiles obtained using the Error Set 2 (Figs 22 and 24) are very similar to those 
obtained with Error Set 1. Obviously, the inversion uncertainties are larger but still very 
similar counter-rotating profiles result. Note that the corresponding y; 2 and AICc values 
(Table]!]) give clear indication of appreciable error over-estimation in the case of Error Set 2 
as already stressed before. Indeed, y 2 < 1 in virtually all cases. This is already evident 
from Fig. [I] where the point-to-point variation of the rotational splittings is visibly much 
smaller than the average uncertainty from Error Set 2. The inversion profiles are equally 
unaffected if the inversion procedure is carried out excluding the most asymmetric splittings 
as evidenced by Figs.[23]and[24j This is not surprising since the large errors associated with 
the most asymmetric splittings already give them less weight compared to the others when 
minimizing y 2 . 


The RLS inversions can be carried out also by regularizing the norm of the first derivative 
of the profile or even the norm of the profile itself, it turns out they are all consistent and 
have similar properties as the RLS inversions using the norm of the second derivative that we 
presented earlier. Figure]!!] shows the corresponding inversion profiles and their associated 
uncertainties derived from the measurement errors. 


So far we have based our inversions on the Model 1 found by Moravveji et al. (2015). 
However, the main qualitative characteristics of the inversion profiles are robust under dif¬ 
ferent model choices. We have explicitly tested that all the models among the best ones from 
the forward modeling in Papics et al. (2014) and Moravveji et al. (2015) produce qualitatively 
similar results for the inverted rotation profiles, i.e., counter-rotation in the radiative enve¬ 
lope. To provide a specific example, we show in Fig. 14 the resulting profile using Model 2 
found by Papics et al. (2014). Although the kernels based on this model do not exhibit a 
large peak right outside the convective core (i.e., no trapped modes), they have a similar 
shape otherwise (see also Fig.]!]). Again, the inversion profile hints at counter-rotation within 
the star’s radiative zone. The uncertainty on the recovered profile is somewhat larger overall 
compared to the uncertainty associated with Model 1 (Fig. 10, top left panel). We attribute 
this to the lack of variation of the kernels in this model, i.e., there are no trapped modes 
to make the kernels more ‘different’ from mode to mode. The more similar kernels for the 
modes of Model 2 result in higher y 2 and AICc values, which we have added for the case of 
RLS and fV = 8 in Table[3] 




















=1 


a-o-4 

iG 


— Norm regularization 

— First derivative regularization 

— Second derivative regularization 


r /R* 


Fig. 13.— The inversion profiles obtained by regularizing the norm of the profile (red, 
/ j -r L s=io-6 , X 2 — 0.53), regularizing the norm of the first derivative (green, /U RL s=io- 6 , X 2 = 
0.51) and regularizing the norm of the second derivative (blue, /i rls=io- 5 , X 2 = 0.52). For 
visual aid, the prohles have been slightly shifted horizontally with respect to each other. 

7. Monte Carlo simulations 

Yet another method to obtain an approximation to the real profile fl(r) under the 
assumption of a smooth profile, consists in generating a large collection of random, synthetic 
test prohles and assigning a score to each measuring how close the predicted splittings are 
to the observed ones. This method is straightforward but very inefficient computationally 
since the number of synthetic prohles that need to be calculated is necessarily large. 

A random, synthetic prohle on a radial grid with Q points can be generated by choosing 
a random rotation value flj at each radial location rj. The rotation values are to be picked 
from a random (uniform) distribution extending from —htoh (in nHz), where the range h is 
to be chosen appropriately as described below. In general, such a prohle will exhibit strong 
huctuations along the radius, i.e. it will be a ’noisy’ spiky prohle, particularly if Q is large 
(we used Q > 100 in our simulations). There are a number of ways to smooth out the prohle, 
a simple one being to use a ‘low pass’ hlter to remove the ‘high frequency’ components of 
the prohle (if we think of it as a time series). The hlter cutoff point defines a characteristic 
length scale A below which the prohle can be considered to have only smooth variations. 
Some padding at each end of the prohle is necessary to avoid end effects when filtering. 
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Fig. 14.— The inversion profile obtained using kernels from Model 2 as found by Papics 
et al. (2014). Here the parameters are identical to those used for the top, left panel in 
Fig. 10 (hrls — 10~ 5 , N = 8). The vertical dashed line marks the convective core radius. 


The counter-rotation within the radiative zone is a characteristic feature of all the models 
we examined. 


See Fig.[l5] for an illustrative example. After the smoothing the profile will have spatial 



Fig. 15.— A random profile (blue) before any smoothing with Q = 100 and h = 500 nHz. 
The ‘low pass’ filtering produces the green curve. The correlation length is A = 0.3. 

fluctuations only on length scales larger than A. Since generally Q > N , i.e., the random 
profiles are defined on a finer grid than the one used to define G, we should resample the 
profile to match the radial grid with N points. We do this by calculating the integral 
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average of the random profile along each segment in the radial grid associated with G after 
an appropriate interpolation. Other methods can be used to achieve the same result. What 
is essential here is that the smoothing should be chosen in accordance to the final resolution 
N so that A N > 1. 

Once the smoothing is performed, the rotation values are not any longer distributed 
uniformly on the [— h, h] interval, resembling instead a Gaussian distribution. This is simply 
because we have introduced short range correlations with our smoothing. Therefore we 
should adjust the range h to ensure that the rotation values at a given radial location are 
more or less equally probable, thus covering uniformly the expected range of fl(r). This 
expected range can be roughly estimated from the mean of the observed rotational splittings 
( 5i/m(3i ). As a concrete example we found that for A = 0.3 and h — 4 /iHz, the random 
profiles visited more or less uniformly the range [—0.2, 0.2] yuHz within a 14% margin. 

If we interpret the rotation value Qj at each radial location as a random variable, and 
given a large collection of random rotation profiles, it is possible to calculate the associated 
covariance matrix. A given row j of this matrix will resemble a Gaussian distribution centered 
at Tj . The mean of the FWHM of the Gaussians in all rows is then an (after-the-fact, of 
course) estimate of A. In practice, the random profiles can be considered approximately 
constant over radial scales not larger than ~ A/3. 



Fig. 16.— Color coded weighted histogram count for each radial location of a large sample 
of random rotation profiles using 1/y 4 as weights (see text for details). The original random 
profiles have been smoothed out so that the typical length scale is A = 0.3. 

Once a random profile has been generated and smoothed out, its associated splittings 
are calculated via Eq. (J2|. We compute then a score which is proportional to x~ 4 (using 
Error Set 1). The lower y 2 is the higher the score becomes. After scoring a large number of 
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profiles (10 9 ), we compute the histogram of the rotation rates at each radial location weighted 
with their corresponding scores and then normalized by an ordinary histogram count. In this 
way, for each radial location and each rotation rate interval we obtain a number indicative 
of its likelihood to explain the observed splittings. 


For the results shown in Fig. 16, the random profiles have a resolution of IV = 14 and 
have been smoothed so that A = 0.3. At each radial segment rj < r < rj + 1 we computed 
the weighted histogram of the ocurrences of Qj (as explained above) over an interval starting 
from -1500 up 1200 nHz and subdivided in 71 bins. By comparing Fig. 16 and Fig. 10 we 
see that the Monte Carlo method reproduces very well the rotation rates at r ~ 0.2 i?* 
while giving only a broad distribution of rotation rates centered around negative values at 
r ~ 0.8 -R*. We note that virtually none of the high scoring random profiles are strictly 
positive (or negative), they all involve at least one sign change along the radial coordinate. 


This use of random profiles is also suitable to establish the ‘quality’ of a set of kernels. 
To do this we take first a random profile (the reference profile) and calculate its associated 
splittings via Eq. (J2]). These splittings are then ‘inverted’ and we compare the resulting 
profile with the reference profile. Some random noise could in principle be added to the 
splittings before attempting the inversion in order to simulate measurement errors, but this 
is unnecessary here since Eq. (14) already describes properly the effect of the measurement 
variance on the inversion profiles. The differences between inverted and reference profiles 
can therefore be attributed solely to the inadequacy of the kernel set to fully recover the 
solution. 
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Fig. 17.— rms error incurred by RLS inversions with Hrls — 10 5 and N = 8 computed 
from random profiles with various degrees of smoothing. See text for details. 


To implement the above, we computed three sets of smooth random profiles (A 









































32 


0.3, 0.4, 0.5, respectively and with 10 6 profiles each). We then rescaled the amplitude of 
each individual profile so as to make the corresponding splittings have a mean that equals 
the mean of the observed splittings. After discarding those profiles whose splittings had 
mixed signs, we proceeded to perform the inversions ( Hrls = 10” 5 ). We computed inverted 
profiles with N = 8 intervals of radial resolution and compared them with the reference pro¬ 
files (integral-averaged over the same radial intervals). The standard errors at each radial 
interval calculated from all the profiles in the set are shown in Fig.[l7j We see clearly that 
the error becomes larger as the profiles have more variability. The radial locations where the 
errors are comparatively smaller coincide roughly with the locations where the A matrices 
have better localization (see Fig. 26). 
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Fig. 18.— RLS inversion profile (solid line) of KIC 10526294 (N = 8, = 10 -5 ). The 

error margins incorporate both the measurements errors (from Error Set 1) and the errors 
related to the inadequacy of the kernels to recover the ‘true’ rotation profile. This latter 
error was estimated by using smooth random profiles with A = 0.3. 


To conclude this Section we present the inversion profile (from the real KIC 10526294 
data) together with an overall (la) uncertainty (derived from both the measurement errors 
and the kernel error as explained above) in Fig. 18 This profile represents a balance between 
good overall statistical measures and good localization properties, at least near the bottom 
of the radiative zone and close to the stellar surface. Note that a fully positive rotation 
profile is possible at 2a level. 
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8. Summary and conclusion 


Numerical models and their pulsation properties (based on the MESA evolution code 
and the GYRE pulsation code) have allowed us to obtain kernels of oscillation modes 
whose frequencies closely match the identified zonal dipole mode frequencies of the B8V 
star KIC 10526294 (Papics et ah 2014; Moravveji et ah 2015). Based on these kernels, we 
computed rotational profiles explaining the detected rotationally split dipole mode frequen¬ 
cies by assuming different functional forms (constant, linear, two-zone and three-zone). We 
also performed RLS and SOLA inversions and implemented a Monte Carlo approach to ob¬ 
tain an approximate rotational prohle and to estimate the errors incurred by the inversion 


process. We relied on the optimal equilibrium model found so far for this pulsator (Moravveji 


et al. 2015) (Model 1) but other seismically derived equilibrium models were also examined 


and lead to qualitatively similar results. 


While the most likely rotation profiles depend on the assumptions made about the 
functional form of the profile, we were able to constrain the average rotation rate near the 
overshoot region to be about 163 ± 89 nHz, a value supported by almost all the rotational 
models we considered. Towards the surface of the star our results are less constrained since 
they are sensitive to the a-priori assumptions on the shape of the rotational profile. If 
a smooth and continuous prohle is assumed, our results point to a mild counter-rotating 
region in the envelope towards the surface of the star rotating at frequency — 717 ± 412 nHz 
with the sign change occurring around r ~ 0.7 i?*. On the other hand, if we assume a 
discontinuous two-zone prohle, we find an outer envelope rotating about six times slower 
than the overshoot region, at 49nHz. The averaging kernel associated with the outer zone 
of this two-zone model leads to a weighted average over most of the radiative envelope. The 
best counter-rotating prohles from inversion, when averaged over the radiative zone using 
this outer zone averaging kernel, lead to rotation rates entirely consistent between the two 
models. 


We performed model comparisons based on the Akaike Information Criterion as well as 
the leave-one-out cross-validation technique, which are both better suited than the reduced 
X 2 when comparing the performance of models that are not-nested, as is the case for the 
models we considered in this study. Both methods give preference to the inversion models 
with the presence of a mild counter-rotation in the radiative envelope at la level. The Monte 
Carlo simulations, fully independent of the above, are consistent with such result. Current 
stellar structure models have so far not considered this type of physical ingredient. 


Following the first rough estimates of H core /f2 e nveiope f° r three core-hydrogen burning 


B stars prior to the asteroseismology space era (Aerts et al. 2003 Pamyatnykh et al. 2004 


Briquet et al.||2007), the recent studies by Kurtz et al. (2014) and Saio et al. (2015) made the 
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first high-precision asteroseismic measurement of surface-to-core rotation in two rv-/ 2.. 5 M 0 
main-sequence hybrid heat-driven pulsators from four years of Kepler photometry. They 
found the star KIC 11145123 to have slightly faster envelope than core rotation and an 
average rotation period ~100d, while KIC 9244992 has slightly faster core than envelope 
rotation and an average rotation period of ~65 d. The authors deduced these results from 
the measured rotationally split g-mode triplets and p-mode triplets and quintuplets without 
relying on forward seismic modeling of the zonal modes as we have done in the present work. 
Our study of the 3.2 M 0 main-sequence B-type star KIC 10526294 hints to an envelope whose 
inner rotation rate is opposite to its outer rate with a small factor ranging from ~ —0.06 to 
~ —0.2 taking into account the uncertainties, while the star has a depth-averaged rotation 
period of about 186 d. In these three cases, even after taking into account that rotation 
rates at stellar birth have been largely overestimated (e.g., Zwintz et al. 2014), a strong and 


efficient mechanism must have been at work to slow down these stars’ rotation after their 
birth. Moreover, an efficient mechanism must be active to transport angular momentum 
within the star. Internal gravity waves (IGWs) could be viable as such a mechanism. Indeed, 
numerical simulations based on IGWs for a 3 M 0 star by Rogers et al. (2013) have shown that 
such waves can transfer angular momentum on short timescales and over the appropriate 
distances in stars with a convective core and a radiative envelope. Additionally, the study 
by Rogers et al. (2013) led to the conclusion that IGWs can lead to either a slightly faster 
envelope than core rotation, or an outer envelope rotating opposite to the inner regions. This 
mechanism thus could be the natural cause of the observational results on the rotational 
properties of KIC 10526294, KIC 11145123, and KIC 9244992. 


The type of rotation profile found for KIC 10526294 and KIC 11145123 is not achieved in 
any standard stellar evolutionary scenario. A similar but much stronger discrepancy between 
models and observations occurs for the core rotation of red giants (e.g., 

In a next step, we do not only plan to perform similar studies as this one for OB-type stars 
with various stellar parameters, but we will also investigate how the stellar structure, and 
in particular the density profile, behaves during the evolution of the star in the presence of 
the most likely rotation profiles we found in this study, testing new physical ingredients such 
as internal gravity waves that were not yet included to describe the physics in the radiative 
envelope of massive stars. Only an extension of the sample of stars with seismic inversion 
treated with appropriate statistical model selection and coupled to an iterative procedure to 
upgrade the input physics can deliver a meaningful improvement in the stellar models. Our 
study is a first step in this direction for massive stars. 


Canticllo et al. 2014). 
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A. Appendix A: Results for Error Set 2 and/or for limited triplet sets 
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Fig. 19.— Same as Fig.[8]with the exception that Error Set 2 has been used. The vertical 
dashed line in the top, right panel marks the location of the minimum of ,\' 2 for the two-zone 
model. 
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Table 4. Same as Table[3]but using the observation uncertainties from Error Set 2. 


Rotation Prohle rms error [nHz] 

V 

x 2 

AICc 

Constant 

11.85 

18.00 

0.45 

12.82 

Linear 

11.78 

17.00 

0.36 

13.67 

Linear+ 

11.66 

16.00 

0.45 

18.05 

Two-zone 

9.61 

16.00 

0.19 

13.93 

Three-zone 

8.44 

15.00 

0.17 

17.13 

RLS, N = 8 

9.17 

15.93 

0.17 

13.83 

RLS, N = 14 

10.29 

16.09 

0.26 

14.66 

SOLA, N = 8 

9.46 

14.30 

0.22 

20.80 

SOLA, N = 14 

7.59 

11.48 

0.15 

35.90 


Table 5. Same as Table[3]but excluding the three most asymmetric splittings, i.e., those 
with periods near 0.96 x 10 5 s, 1.29 x 10 5 s, and 1.83 x 10 5 s. The single m = +1 mode with 
period near 1.78 x 10 5 s was also excluded. Error Set 1 is used. 


Rotation Prohle 

rms error [nHz] 

V 

x 2 

AICc 

Constant 

9.07 

14.00 

25.07 

355.9 

Linear 

7.22 

13.00 

8.55 

119.4 

Linear-I- 

8.29 

12.00 

22.08 

277.0 

Two-zone 

4.95 

12.00 

1.69 

32.27 

Three-zone 

5.16 

11.00 

0.53 

22.46 

RLS, N = 8 

4.84 

11.29 

0.63 

22.33 

RLS, N = 14 

4.77 

11.20 

2.35 

42.00 

SOLA, N = 8 

5.86 

9.45 

0.86 

34.48 

SOLA, N = U 

3.58 

7.21 

0.40 

53.45 
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Fig. 20.— Same as Fig.[19]with the exception that the three most asymmetric splittings have 
been excluded from the y 2 minimization, i.e. those with periods near 0.96 x 10 5 s, 1.29 x 10 5 s 
and 1.83 x 10 5 s. The single m = +1 splitting with period near 1.78 x 10 5 s was also excluded. 
Error Set 1 has been used. 
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Fig. 21.— Same as Fig. 20. i.e., omitting the three most asymmetric modes and the single 
m = +1 splitting, but this time using Error Set 2. 

Table 6. Same as Table[5]but using the observation uncertainties from Error Set 2. 


Rotation Profile rms error [nHz] 

V 

X 2 

AICc 

Constant 

8.95 

14.00 

0.47 

11.63 

Linear 

7.05 

13.00 

0.32 

12.36 

Linear+ 

8.27 

12.00 

0.47 

17.66 

Two-zone 

4.91 

12.00 

0.14 

13.67 

Three-zone 

4.70 

11.00 

0.13 

18.11 

RLS, N = 8 

4.96 

12.15 

0.14 

13.11 

RLS, N = 14 

5.85 

12.39 

0.24 

13.42 

SOLA, N = 8 

6.04 

10.99 

0.23 

19.27 

SOLA, N = 14 

3.49 

8.59 

0.10 

34.55 
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Fig. 23.— Same as Fig. [IQ] but this time excluding the three most asymmetric splittings, i.e., 
those with periods near 0.96 x 10 5 s, 1.29 x 10 5 s and 1.83 x 10 5 s. The single m = +1 mode 
with period near 1.78 x 10 5 s was also excluded from the computation. 
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Fig. 24.— Same as Fig. 23 with the exception that Error Set 2 has been used. 
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Fig. 25.— Same as Fig. 11 but excluding the three most asymmetric splittings, i.e., those 
with periods near 0.96 x 10 5 s, 1.29 x 10 5 s and 1.83 x 10 5 s. The single m = +1 mode with 
period near 1.78 x 10 5 s was also excluded from the computation. 
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B. Appendix B: The A matrices 


As mentioned in Section [6j the matrix A gives an indication of how well the inversion 
profile recovers the true profile (in the ideal case of no measurement error in the splittings). 
Using the best model from Moravveji et ah (2015), we can see this fact at work very clearly 
in both Figs. 26 and 27 as the respective fi parameter is varied. The more A resembles the 
identity matrix the better the reconstruction is, thus providing a qualitative assessment of 
the inversion. 
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Fig. 26.— The matrix A for the RLS method using the kernels from the best model in 

with two different resolutions (N = {8,14}; top and bottom row 
respectively) and a range of smoothing parameters Hrls■ Each row on these matrices is the 
discrete version of the averaging kernel Jf(r',r). The A matrix should resemble as much 
as possible the identity matrix, as Eq. ( [I3| ) indicates, if we want a faithful reconstruction of 
the rotation rate (at the radial location defined by the row index, i.e., first row corresponds 
to the first bin in the radial grid). Note how near r/R* ~ 0.17 the rotation rate is well 
recovered almost independently of the choice of N or /irls- 


Moravveji et al. (2015 
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VSOLA 
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Fig. 27.— Same as in Fig. 26 but using the SOLA method. The matrices are very well 
localized except at higher resolutions and higher /j,sola • At / isola = 10 -2 and N — 8 there 
is a small amount of leakage from regions close to r = 0.4/?* into the outer radial bins 
(r' ~ 0.9/?*). We estimate from this that the rotation rate there is overestimated by the 
inversion by about 3%. 


C. Appendix C: Averaging kernels in the continuous limit 

The RLS and SOLA inversions behave very differently when the resolution is increased. 
In this section’s experiments we kept /irls = ICO 5 and [isola = 10~ 2 fixed but used three 
different resolutions, N = 14,28,5000. From the viewpoint of the inversion, it is possible 
to increase the resolution beyond the number of observations M since regularization keeps 
the effective number of fitted parameters below M. Here we do not include the variances on 
the inversions in the discussion. We have chosen the Error Set 1 as the uncertainties on the 
splittings. 

In Fig.[28]we show the resulting inversion profiles. There are two radial locations where 
most of the rotation rates roughly coincide. One is at r ~ 0.17 and the other one is at 
r ~ 0.92 where the rotation values are not too far from each other except for the SOLA 
inversions with N = 28, 5000. 

Let us examine the averaging kernels from three selected radial locations r 0 = 0.17 /?*, 
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Fig. 28.— Rotation profiles from RLS (^rls — 10 -5 ) and SOLA ( Usola = 10 -2 ) methods 
and three different resolutions N = 14, 28, 5000. The inversions are based on the kernels 
from the best model of Moravveji et al. (2015). The measurement uncertainties are taken 
from the error set 2. 


r 0 = 0.48 R* and r = 0.92 R* as shown in Figures [29j [30] and [31J respectively. From the figures 
we see that for r 0 = 0.17R* all the averaging kernels are indeed well behaved generally, so we 
expect inferences for this location to be consistent. At r 0 = 0.48 (Figure |30|) the localization 
is acceptable as long as N is low. Closer to the stellar surface, at r 0 = 0.92 (Fig. 31), the 
situation is similar although the RLS kernels degrade considerably already when N = 28. 


From these figures we conclude that consistent inferences are to be found using either 
RLS or SOLA methods if the resolution is kept low, i.e., N < 14. Such a resolution also 
implies a sensible model comparison through the AICc (see Table[3]), while larger N and small 
smoothing parameters would lead to an effective number of degrees of freedom v < 1 and 
strongly negative AICc values, implying overfitting from the viewpoint of model comparison. 
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Fig. 29.— The row of A corresponding to r 0 = 0.17 at three resolutions (top three plots). 
Bottom plot is a zoom of the one immediately above. 
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Fig. 30.— Same as Figure [29] but for r 0 = 0.48. 
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Fig. 31.— Same as Figure [30] but for r 0 = 0.92. 
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D. Appendix D: Testing RLS inversions with a synthetic profile 

The following test is to check that the counter-rotation profiles are not produced by some 
undesired property of the RLS inversion methods. We take the optimum two-zone model 
from Section 5 and smooth it using the ‘low pass’ filter with correlation length A = 0.3 
described in Section 7. This is taken as the ‘actual’ rotational profile, which does not exhibit 
counter-rotation. Subsequently, we calculate the associated exact rotational splittings via 
Eq.@. To each of these 19 splittings we add random noise sampled from a Gaussian 
distribution with zero mean and the same standard deviation as the actual measurement 
errors (Error Set 1). We set N = 8, ls — 10 -5 as used for the RLS inversions in the main 
text and proceed to calculate the inversion profile. 

At each radial bin we compare the inversion value with the integral average of the 
‘actual’ profile over the same radial bin. This gives us a direct estimate of the inversion 
error. By repeating this process a large number of times we can obtain well defined statistics 
(we used 10' iterations). Figure[32| shows the ‘actual’ profile in blue, the recovered profile in 
black and the estimated la uncertainty range in red. The RLS method does a good job in 
recovering the actual profile, which always occurs within the errors. Some inversion profiles 
must counter-rotate mildly since the error region extends below zero in the outer half of the 
star. However, given the errors, we do not find a counter-rotating profile in this case. 



Fig. 32.— The smoothed-out two-zone model (blue, correlation length A = 0.3) and the 
recovered profile from RLS inversion (black, /trls = 10~ 5 ). Error bounds are in red. 
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